#ifndef PARTIALDER_H
#define PARTIALDER_H

#include "Core/RowSparse.h"
#include "Embedded/EmbeddedGrid.h"

template <int Dim, int pD, int _fitOrder = 5>
class PartialDer
{
public:
  using iVec = Vec<int, Dim>;
  using rVec = Vec<Real, Dim>;
  using EmGrid = EmbeddedGrid<Dim>;
  enum { fitOrder = _fitOrder };

  struct partialDerTransform;

  PartialDer(const EmGrid &aGrid);

  void apply(const Tensor<Real, Dim> &aData, Tensor<Real, Dim> &outData) const;

protected:
  const EmGrid    &emGrid;
  RowSparse<iVec> matOp;
  const Real      dangerousCN = 1e4; // give a level-1 warning if the condition number is larger than this.
};

//============================================================

template <int Dim, int pD, int _fitOrder>
struct PartialDer<Dim, pD, _fitOrder>::partialDerTransform
{
  bool checkStandardScheme(const Tensor<int, Dim> &nodeLabels, Vec<int, Dim> idx, bool useOuterRect) const
  {
    const int radius = (fitOrder-1) / 2;
    auto validBox = nodeLabels.box();
    if(idx[pD] < validBox.lo()[pD] + radius || idx[pD] > validBox.hi()[pD] - radius)
      return useOuterRect;
    auto u = Vec<int, Dim>::unit(pD);
//    using NodeType = typename EmbeddedGrid<Dim>::NodeType;
    for(int i = -radius; i <= radius; ++i)
      if(nodeLabels(idx + u*i) < 0)
        return false;
    return true;
  }
  template <class... Ts>
  Real transform(const Vec<Real, Dim> &p, Ts... indices) const {
    Vec<int, Dim> beta { decltype(indices)::value ... };
    if(beta[pD] == 0)
      return 0.0;
    --beta[pD];
    Real r = 1.0;
    for(int d = 0; d < Dim; ++d)
      r *= std::pow(p[d], beta[d]);
    return r * (beta[pD]+1);
  }
};

#endif //PARTIALDER_H
